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Abstract 

The red blood cells or erythrocytes are biconcave shaped cells and consist mostly in a membrane 
delimiting a cytosol with a high concentration in hemoglobin. This membrane is highly deformable and 
allows the cells to go through narrow passages like the capillaries which diameters can be much smaller 
than red blood cells one. They carry oxygen thanks to hemoglobin, a complex molecule that have very 
high affinity for oxygen. The capacity of erythrocytes to load and unload oxygen is thus a determinant 
factor in their efficacy. In this paper, we will focus on the pulmonary capillary where red blood cells 
capture oxygen. We propose a camera method in order to numerically study the behavior of the red 
blood cell along a whole capillary. Our goal is to understand how erythrocytes geometrical changes 
along the capillary can affect its capacity to capture oxygen. 

The first part of this document presents the model chosen for the red blood cells along with the 
numerical method used to determine and follow their shapes along the capillary. The membrane of 
the red blood cell is complex and has been modelled by an hyper-elastic approach coming from 1 16 1. 
This camera method is then validated and confronted with a standard ALE method. Some geometrical 
properties of the red blood cells observed in our simulations are then studied and discussed. 

The second part of this paper deals with the modeling of oxygen and hemoglobin chemistry in the 
geometries obtained in the first part. We have implemented a full complex hemoglobin behavior with 
allosteric states inspired from |4|. 

Les globules rouges ou erythrocytes sont des cellules biconcaves et sont formees d'une membrane 
contenant le cytosol, fluide ayant une grande concentration en hemoglobine. Cette membrane est tres 
deformable et permet a ces cellules de traverser des passages etroits comme les capillaires dont les 
diametres sont souvent plus petits que ne Test le diametre des globules rouges. Les erythrocytes trans- 
portent l'oxygene dans 1'organisme grace a la proteine d'hemoglobine qui dispose d'une grande affinite 
pour l'oxygene. La capacite des globules rouges a capturer et liberer l'oxygene est ainsi un facteur 
determinant de leur efficacite. Dans cet article, nous nous restreindrons a l'etude d'un capillaire pul- 
monale ou les globules rouges capturent l'oxygene. Nous developpons une methode de camera de facon 
a etudier numeriquement le comportement de globules rouges le long d'un capillaire. Notre but est de 
mieux comprendre comment les eythrocytes se deforment au long du capillaire et comment cela affecte 
leur capacite a capturer l'oxygene. 

La premiere partie de ce document presente le modele que nous avons choisi pour representer le 
globule rouge, ainsi que la methode numerique utilisee pour determiner et suivre la forme du globule le 
long du capillaire. La membrane du globule rouge est complexe et a ete modelisee par une approche 
hyper-elastique issue de 1 16 1. Cette methode de camera est ensuite validee et confronted a une methode 
ALE standard. Enfin, quelques proprietes geometriques du globule rouge resultant de nos simulations 
sont etudiees. 

La deuxieme partie de ce papier traite de la modelisation de la chimie entre l'hemoglobine et l'oxygene 
dans des geometries issues de la premiere partie. Un modele complet du comportement de l'hemoglobine 
a ete implements, en particulier faisant intervenir deux etats allosterique des molecules d'hemoglobine. 
Ce modele s'inspire de |4). 
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Introduction 



Human respiration and more generally mammal respiration, is a complicated efficient system that uses 
oxygen amongst other products in order to provide energy to the body. One of the function of the respiration 
is thus to bring oxygen from air into the mitochondria where the major energy source of the body, the ATP, 
is synthesized. Oxygen travels a complex exchange system in different media. Ambient air is first brought 
on an air/blood exchange surface (the aveoli) through the lung, where oxygen is extracted from air. The 
fractal like structure of lungs (dichotomous tree of twenty-three generations) permits the exchange surface 
to be large enough (~ 100 to 2 ) to fit energy needs of the body. The particular structure and functioning 
of the lung are very important in respiration and have been intensively studied, see for example (9][15 | but 
mostly independently of blood behavior. The next step of respiration for oxygen is to cross the exchange 
surface, to reach the red blood cells in the capillaries (smaller blood vessels) and to react with hemoglobin 
in order to be carried away from the lung towards the cells and mitochondria through an other branched 
complex structure, the vascular network, which also has been studied intensively, see for example ifTUll . 

In this paper, we are interested in the modeling of the step of oxygen transport that takes place in 
the pulmonary capillaries. There, oxygen crosses the membrane between the alveola and the capillary and 
dissolves into plasma, where oxygen is transported by a convection/diffusion process towards the red blood 
cells. Next, oxygen crosses the membrane of the red blood cell and diffuses inside it (in the cytosol), then 
meeting an hemoglobin molecule, a chemical reaction can occur and oxygen is trapped by the hemoglobin 
EJ. 

Such a modeling can lead to a better understanding of the phenomena involved in oxygen exchanges at 
capillary level. In particular, some pathologies can affect these exchanges by modifying critical parameters. 
For example, the rigidity of the red blood cells membrane can be modified (paludism, diabetes) and lead 
to less efficient shapes in terms of capacity to cross capillaries or capacity to catch oxygen. Similarly, the 
chemical properties of hemoglobin are very sensitive to blood parameters (as pH) and global consequences 
on oxygen transport in blood could be determined. 

The convection/diffusion process involved in oxygen transport in the capillaries is highly dependant on 
the shape of the structures (capillary and red blood cells) but also on the velocity field of the plasma. Hence 
the determination of these are required in order to correctly model oxygen displacements. The capillaries 
diameters are often much smaller than the red blood cells' ones, but thanks to their highly deformable 
membrane the red blood cells shapes can change under external constraints and they can cross the smaller 
vessels (parachute shape). 

In order to deal with these different points, a numerical camera method is used to follow the red blood 
cells deformation along a straight capillary. This method is based on the standard ALE method for fluid 
structure interaction but is written in a moving frame. The benefits of this approach is that it is not limited 
by mesh skewing and the whole capillary can be covered. The mechanics of the membrane of the red blood 
cells is modeled by a Yeoh hyper-elastic model from [16] and validated. Some properties of red blood 
cells (velocities, shapes in the capillaries) are studied along with apparent hydrodynamical resistances of 
capillaries crossed by red blood cells. 

The second step of this paper is to add to the previous modeling the behavior of oxygen in the plasma 
and in the cytosol, and the behavior of hemoglobin in the cytosol. Two physical processes are involved: 
displacement (convection and diffusion for oxygen, diffusion only for hemoglobin) and chemistry (oxygen 
and hemoglobin reactions). A realistic model of hemoglobin chemistry is used [4| and implemented nu- 
merically. Some preliminary results are presented, mostly on the consequences of capillaries diameter on 
the red blood cells capacity to capture oxygen. 

Note that all simulations have been performed using Comsol Multiphysics 3.3a and that most of our 
calculations are 2D-axisymmetric. 
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1 Fluid structure interaction: modeling 



1.1 Geometries 
1.1.1 Red blood cells 

The red blood cell can be considered as vesicle with a membrane enclosing a fluid, the cytosol. The shape 
of the red blood cell is particular and known to minimize the bending energy of the membrane |2| amongst 
other properties ll22l . see figureQ] In the following work, the geometrical parameters of the red blood cells 
have been chosen accordingly to Weibel's data ll23ll for human erythrocytes. 

The cytosol can be considered as a viscous fluid and contains a high concentration of hemoglobin. It 
will be considered as a Newtonian viscous fluid: its density corresponds to that of water while its viscosity 
is eight times water viscosity [|T4l . 

The membrane of the red blood cells is a superposition of a bilipidic layer and of an actin mesh. Its 
thickness is of the order of one or two hundreds of nanometers |[23ll . 




Figure 1: Geometrical model used for the red blood cells. The radius of the red blood cell is 7.34 fim, its 
thickness ranges from 1.4 to 2.4 /im (data from Weibel ]23l ). The membrane thickness is 200 nm. 

The mechanical behavior of the membrane is complex and not yet well understood. However some 
interesting models have been developed |fl] [16] QjD . More precisely, in fl6l . a relatively simple hyper- 
elastic model for the membrane has been applied successfully to recreate numerically optical tweezers 
experiments. This hyper-elastic model has been developed by Yeoh in [25 | and uses the following strain 
energy: 

W = ^f(Xi + \l + X 2 3 -3) + C{\\ + + 3) 3 (1) 

where Go is the initial shear modulus of the membrane and G a constant estimated to Go/30 in order 
to fit optical tweezers experiments, see ifTBI . Note that in the literature, data are often given in term of 
in-plane shear modulus: /io = Goho where ho is the initial thickness of the membrane. In the following, 
/io has been chosen to 11.3 fiN/m, see ll.1.21 

1.1.2 Reproducing Mills et al experiments 

In order to validate the membrane model, a numerical reproduction of Mills et al experiments [ 16 1 has been 
performed. These experiments consist in attaching a metallic bead at each extremity of the red blood cell 
and to use magnetic sinks in order to control the positions of the beads (optical tweezers). This method 
gives access to the force applied to the beads and thus it is possible to deduce some mechanical properties 
of the membrane if measuring the deformations against the force applied. 

In order to build a numerical equivalent, we use our three dimensional red blood cell. Two cubic rigid 
"beads" have been attached on its membrane. One of these beads will be fixed while the other will be 
submitted to a force. The inside of the red blood cell is considered as a Stokes fluid with previously given 
characteristics. Note that because the fluid is assumed incompressible, the red blood cell keeps a constant 
volume. We use an ALE fluid-structure interaction model. Figure[2]shows the initial geometry and a typical 
deformed globule under stress. 

We have measured lengths, widths and thicknesses of the red blood cell for a range of forces going 
from to 210piV with steps of 30pN. We have found a better collusion between our results and Mills 
et al results for an in-plane shear modulus of 11.3/j.N/m and with the parameter G equal to Gq/30. The 
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Figure 2: Left: initial globule geometry, with "beads" at each side. Right: deformed globule with a force 
of 120pN. 




Force (pN) 



Figure 3: Deformation of the red blood under optical tweezers stress, numerical results. 



corresponding results of our numerical simulations are shown on figure [3] In the following, we will keep 
these parameters constant to these validated values. 

1.1.3 Capillary 

We will assume at first approximation that a capillary is topologically equivalent to a cylinder whose axis 
can be considered as a straight line. In the camera method developed below, this will simplify the movement 
of the camera, limiting it to a translation along the capillary axis. The walls of the capillary will be assumed 
unmoving and rigid. The shape of the capillary will be chosen according to measures done in EH . The 
typical length of such a vessel is of the order of one millimeter and its radius can range from h[im to 15/Ltm. 
Figure|4]shows a geometrical model of capillary that will be used in the next sections. 

The plasma will be considered as a viscous Newtonian fluid with same density than water, but with 
twice its viscosity. Because capillaries are far from the heart, they do not feel heart beating and thus 
we can consider that the plasma circulates with a constant flow rate along the capillary. Moreover the 
velocities are of the order of the millimeter per second and the Reynolds number pvdj '/i can be estimated 
to 1000 x 10~ 3 x 10~ 5 /2.10~ 3 — 0.005. Thus Stokes equations are a good approximation of plasma 
behavior in the capillary. 



Figure 4: Example of a capillary 
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1.2 Equations 



Because red blood cells are travelling through the whole capillary, the coordinates map induced by their 
displacement in a standard ALE method will be highly skewed. Direct numerical simulations using finite 
elements will thus be very difficult because of highly deformed elements. This could be bypassed by 
regular remeshings but this method is relatively costly and relatively difficult to set up. 

Hence, considering the typical shapes of the studied capillary (that is assumed to be topologically 
equivalent to a cylinder and to keep its axis along a straight line) and considering the fact that our objects 
of interest are the red blood cells, a camera method seems advantageous. This method consists in focusing 
only on the portion of the capillary where the red blood cells are and to displace focus together with red 
blood cells. 

Numerically speaking, this consists in meshing a small portion of the capillary, in making the mesh 
move along it at red blood cells mean velocity and also in making the mesh fit capillary local shape. Note 
that because we do not need to mesh the complete capillary, this method will also benefit of a reduced 
number of elements. 

The drawbacks of this method are located on boundaries: 

• the coordinates map induced by red blood cells deformations must also include geometrical changes 
of the capillary walls along its length 

• the boundary conditions for the plasma on the non-walled boundaries are not obvious and need to be 
studied more closely 

Note that all bearings involved in the following will be built thanks to the Winslow smoothing method 
l24l . which robustness against mesh elements inversion is better than Laplace method. 

1.2.1 Standard method 

First we will write the standard method, namely the standard fluid structure equations in ALE formulation 
on the full capillary. The capillary domain will be called O and its straight axis will correspond to one of 
the reference frame axis. In the following, this particular direction will be recognized amongst the other by 
the subscript c. The reference coordinate system will be denoted (x), and the basis vectors will be noted 
(ej) j. Thus the cylinder axis is directed by e c and the position on this axis referenced by x c . 

Let G C £1 be the membrane of the red blood cells and F = Q\G the reference domain for the fluid. 
Note that F corresponds not only to the outside of the red blood cells (plasma), but also to their insides 
(the cytosol), these two fluids have the same density (similar to water) but a different viscosity. Let 7$ 
(i = 1, 2, 3) be the capillary walls as shown on figure (0 such that d£l — 71 U 72 U 73. 



G 




Yi 





Figure 5: Reference geometry for standard ALE method. 



The standard fluid structure equations in ALE formulation are: 



p g y-div(a(u))=0 



onG 




(2) 
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with boundary conditions: 



cr(u).n= a ip 1 {nfS7v.n — p.n)" ondG 

v =dg ondG t = i,{dG,t) 

v = on ^(71, f) = 71 

v = v e on -0(72, t) = 72 

/z/Vw.n — pn = on -0(73, =73 

and ip = Id + p, where p is a bearing of u on whole fi with boundary constraints i/p = u(., t) on <9G 
and p.n = on 3S1. The constraints equality is done on the reference frame and the notation "ip~ x (...)" 
means that the constraints tensor of the fluid defined in the deformed frame is pulled back in the reference 
frame. The boundary condition on dG is the displacement due to the deformation of the membranes of the 
red blood cells, this deformation comes from equations (0. The boundary condition on dfl states that the 
capillary walls are unmovable but allows the coordinates of the deformed frame to slide along the walls of 
the capillary. The coordinate in the deformed frame ip(fl, t) will be noted y and y — ip(x, t). 



1.2.2 the camera method 

We will now restrict the domain for the camera method: we will now consider a coordinates system (X) 
directed by (Ei)i and assume our reference frame is a cylinder C which axis is along X c . Its length along 
X c is L > (L is smaller than capillary length). This cylinder will fit the local capillary geometry while it 
follows the red blood cells along the capillary. The deformed cylinder at time t will be called the camera 
and will be noted C(. The domain defined by the red blood cells membranes is called G C C. 

The camera will follow the red blood cells in such a way their mean position is always at the center of 
the view of the camera. Hence the position along the capillary of the camera center can be calculated: 

mc{t)=m o + kMA (4) 
Jg 1 

The displacements u are still calculated in the reference frame on the membrane domain G C C 
and will not differ from the standard ALE method. In order to rebuild the geometry of the capillary 
containing the deformed globules, it is necessary to first integrate the deformed globules into the non- 
deformed capillary and then to make the cylinder wall fit the capillary without changing the shapes of 
the globules. This will be done using two successive bearings Y = A(X,t) = X + a(X,t) and Z = 
B(Y,t) = Y + (3(Y,t): 




Figure 6: Domains notations and bearings involved in the camera method. 



Thanks to the definition of m c (t), the transformation A preserves the mean position of the red blood 
cells to on the X c axis. 



6 



Stokes equations (velocity will be called w and pressure p) will be calculated in (Z) coordinates in order 
to be solved in the correct geometry. The transformation A corresponds to the movement of the camera 
and will thus induce an additive transport term in the Stokes equations: — (^jp.V)u. However, because of 
the small Reynolds numbers involved and thanks to the fact that ~ v, this supplementary term is 

negligible the same way (v.V)v is. Moreover, the transformation B does not modify Stokes equations. 

Finally, the succession of these two bearings can be reduced to a single bearing S(X) = X +£(X, t) — 
X + a(X, t) + (3(X + a(X, t),t). The boundary conditions for £ are: 

£ = «(.,*)- m c (t).E c ondGt 

Cn(X)=g(X + aX,t) + m c (t).E c ) on T 1 (6) 
£.n = on T 2 and T 3 

The boundary conditions for fluid equations will be the same as for the standard ALE method con- 
cerning the capillary wall and the globules wall (non slip conditions). However, they should be studied 
for what concerns the inlet (H(ra)) and the outlet (S(T3)) of the camera domain. In the standard ALE 
method, a Dirichlet condition on velocity was imposed at the outlet 73 and a zero pressure condition was 
introduced at the inlet 72. The Dirichlet condition can be easily used at camera outlet considering the 
very low Reynolds numbers: a Poiseuille profile should be a good approximation of the velocity at camera 
outlet, all the most if the border is not too close to the red blood cells. At inlet (Z c = —L/2) however, 
the correct pressure condition should be p(Z, t) — p(y, t) with Z such that Z c = —L/2 and y such that 
y c = m c (i) — L/2. Again, thanks to the low Reynolds numbers, we can assume that pressure is constant 
on the section y c = m c (t ) — L/2 and impose p to fit the mean value on this section namely po(t). However 
this value is only time dependant and because Stokes equations only depends on Vp, making the variable 
change q(Z, t) — p(Z, t) — po(t) will not change the equations and will allow us to impose q = on the 
camera inlet. 

Finally, the fluid constraints on red blood cells membranes — p.n + fifWu.n is modified only by the 
relative pressure change: because this change is global it does not affect the mechanics equations. 

The errors done relatively to the standard approach are due to boundary conditions simplifications and 
to the suppression of the new transport term. We will do numerical comparisons in order to show the quasi 
equivalence between both methods and show that, in the contrary of the standard approach, the camera 
method can follow the red blood cells all along the whole capillary. 

2 Fluid structure interaction: numerical results 

2.1 Comparison between the standard method and the camera method (2D ax- 
isymmetric) 

In order to validate the camera method, a comparison has been made between the standard ALE method 
and its pendant in camera method. In order to limit the calculation time, a 2D axisymmetric case has been 
used. In this example, a red blood cell has been put inside a plasma flow through an axisymmetric bottle- 
neck. At first the radius of the capillary is larger (4.5/xm) than that of the globule (3.6/im) but it is smaller 
in the constriction (3/im), see figure[7] 



Figure 7: Capillary geometry used in this section, the red blood cell enters to the left part where radius is 
4.5pm, the full structure is 1.5mm long and the constricted part is 0.5mm long with a radius of 3fim. 

Because of its design, the camera method induces a far smaller distortion of the mesh and can follow 
easily the red blood cell along the whole geometry. The mesh quality for the camera method is only de- 
pendant on the radius of the capillary and not on the position of the red blood cell. On the contrary, mesh 
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quality is decreasing linearly in the standard method, and the numerical simulation stops before the end of 
the bottleneck (at time 1.74s about 600 [im from entry) because of highly skewed elements. Mean mesh 
qualities versus globule positions are represented on figure [8] 



Camera 
Fixed 



500 1000 
Red blood cell displacement (micrometers) 



Figure 8: Comparison between standard ALE method (red dashed curve) and camera ALE method (black 
continuous curve) in a bottleneck like "capillary": mesh qualities. 

Some parameters have been checked in order to compare the differences between both methods. We 
have measured the displacement of the red blood cell and also its aspect ratio along time in both simula- 
tions. The data show a very good agreement between both methods on the range of convergence of the 
standard method. There is however a slight difference for the red blood cell aspect ratio in the last times 
of convergence of the standard method that is a consequence of some of its elements being highly skewed: 
this lead to a bad convergence of the problem solution. 





Figure 9: Comparison between standard ALE method (red dashed curve) and camera ALE method (black 
continuous curve) in a bottleneck like "capillary". Left: position of the red blood cell in the capillary. 
Right: length over diameter ratio of red blood cell. 



2.2 Going along a capillary 

Being able to follow red blood cells along a capillary gives access to a large quantity of data that needs to 
be studied in lots of different contexts. We will present here some of these interesting results but keeping 
in the wide context: "how do it feel being a red blood cell in a capillary ?". 

In this section, and as before, we will limit us to a 2D axisymmetric problem in order to limit calcula- 
tions times. 
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We will do our study on an infinite capillary along axis x c (measures will be done in meters). The 
vessel diameter is equal to 8.5^m on x c £ ] — oo, 0[ U ]170.10 -6 , +oc[. and is variable on the range 
x c 6 [0, 170. 10~ 6 ]. The variations chosen are close to the one observed on an in-vivo rat capillary in 
[1 1 31 and are shown on figure [TT)1 We will call C, the nortion of canillarv corresnonding to the region of 
diameters variation, nam« 

5 1 | i | i | i | i | i | i 
4.5 - 



1 4 - 
i 3.5 - 




-50 50 100 150 200 250 

Position in capillary (micrometers) 

Figure 10: Capillary radius. 

Depending on the case, we will put one, three or five red blood cells in the capillary. The initial distance 
between them has been put to 10[im in order to avoid collisions. This value is higher than the observations 
that gives a number closer to 5 [im (data calculated from a red blood cells volume fraction of 0.4 in blood, 
see |fl9l 1231 ). So that they have a stabilized shape at the entry of the portion of capillary C v , the initial 
mean position of the globules has been put to x c = — 650/xto. Plasma flow through a capillary section of 
8.5/xm diameter has been chosen to correspond to a mean velocity of 0.5mm/ s, in the simulations the flow 
is maintained constant whatever the position of the camera, in particular it implies velocity increase when 
capillary radius decreases. 

2.2.1 What shape? 

Under plasma flow stress, the red blood cells are moving forward the capillary, and because of the non 
homogeneous shape of the flow on a capillary section, the globules shape changes to the typical parachute 
shape, pointing in the direction of the flow. Moreover the globules go through skewer shapes when the 
radius of the capillary decreases. In order to estimate the distortion of red blood cells, we will use the 
aspect ratio as in |fl3l . It is defined by the ratio of cells length along the capillary axis over their diameter, 
see figure [TT1 Under no stress, this aspect ratio is 2.4/7.34 ~ 0.33. 



Figure[TT]shows the different curves of aspect ratio depending on the number of red blood cells involved 
in the calculation. Red blood cells aspect ratio follows radius changes of capillary. This illustrate how able 
to deform the red blood cells are: starting from an aspect ratio of 0.33, they can reach almost 2 in the 
narrowest parts of the capillary. These numbers are fully coherent with the measures made in (T3), where 
the authors have observed similarly the red blood cells aspect ratios and have more particularly studied 
these changes against capillary radius variation. On figure [T2l typical shapes observed in our simulations 
have been drawn, in particular we can observe the classical parachute shape. 

The presence of more than one red blood cells plays a small but visible role in aspect ratio determination 
for our choice of parameters (distance of 10/im between each cell). This indicates that closer globules, as 
in real blood, should influence themselves, all the more if there are collisions between them. 

2.2.2 Which velocity? 

Because the skewed red blood cell does not fill the whole section of the capillary, it is submitted only to 
the larger part of plasma velocities which stands mostly in the middle of the capillary (Poiseuille profile). 
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L Position in capillary (micrometers) 



Figure 11: Length (L) over diameter (D) ratio of the middle red blood cell. 




Figure 12: (a) Initial shape of the red blood cell. (b),(c),(d) Shapes of red blood cell under flow stress for 
different values of capillary radius. 



Thus, the mean velocity of the cells is most of the time larger to the mean velocity of the plasma. This 
characteristic of blood is well known [7]. Similarly, it is interesting to observe that the gap between the 
cells and the capillary wall decreases with capillary radius and can indeed modify this phenomena. The 
velocities of plasma and the velocities of cells in different cases (one, three or five cells in the capillary) are 
shown on figureQj] Note that, for the left figure, increasing the number of cells in the capillary will induce 
an homogenisation of the mean velocity because each cell sees a different radius than the others. On the 
right part of the figure the velocities of the middle red blood cell is shown, the effects of the presence of 
other cells is small but exists. Hence, decreasing the initial distance between cells should induce a higher 
interaction, but also create collisions which are not, for now, managed by our code. 



2.2.3 Consequences on the hydrodynamical resistance of the capillary? 

By crossing the capillary, the red blood cells interfere with vessel's hydrodynamical resistance which can- 
not be predicted through fluid models, even non-newtonian ones. This is due to the deformability of the 
red blood cells along with the similar sizes of red blood cells and capillary diameters. Thus it is necessary 
to study the full fluid structure interaction problem in order to obtain reasonable high scales estimation of 
capillaries hydrodynamical resistances. In order to apply our simulations at this problem, we have obtained 
preliminary results on the hydrodynamical resistance. We have measured the resistance of the portion of 
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Figure 13: Left: mean velocities of the red blood cells. Right: velocity of the middle red blood cells. 

capillary viewed by the camera against the number of cells (zero, one, three or five). To estimate resistance 
we have calculated the ratio from pressure drop between the inlet and the outlet over the flow in a capillary 
section. 

The results are drawn on figure[l4]and show that, as expected, increasing the number of cells increases 
the resistance. Note that at relative high diameters, the resistance is almost identical whatever the number 
of the globules. On the contrary, the differences increases when globules shape is changing, actually some 
of the fluid energy is exchanged with the membrane of the red blood cells and higher pressure drop is thus 
needed to ensure constant flow 



1 red blood cell 
3 red blood cells 
5 red blood cells 
no red blood cell 




Position in capillary (micrometers) 



Figure 14: Resistance of the portion of capillary in the camera versus camera position. 



3 Introducing oxygen and hemoglobin: modeling (preliminary) 

Red blood cells are designed to catch oxygen in the lung capillaries, to carry it in the whole body and to 
deliver it to the cells in the systemic capillaries. The way oxygen is captured and delivered is obviously 
critical. Actually, the time spent in the capillaries must be tuned in order to maintain the balance between 
oxygen delivery and oxygen needs. The chemistry of hemoglobin plays an important role in oxygen ex- 
changes, but other actors are also involved, like blood velocity or red blood cells shapes that are determinant 
for good functioning of the system. 

In this section we will use the QD hemoglobin chemistry model developed in [4| and extend it to 2D 
axisymmetric one in the geometries obtained from the preceding sections. More precisely and because 
oxygen diffusion and reaction with hemoglobin can be uncoupled from red blood cells deformation, the 
results of the preceding sections can be considered as known data in the modeling of oxygen exchanges. 
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3.1 Modeling of hemoglobin chemistry 



Hemoglobin is a complex molecule consisting in four hems, each of them being able to carry one molecule 
of oxygen. The four hems divides in two a and two (3 hems (or chains) that have different bindings 
properties with oxygen. This asymmetry of the hemoglobin molecule implies the existence of two allosteric 
states T and R which will be integrated in the model. It is not clear in which priority these chains react 
with oxygen, and we will follow one of the mechanisms proposed in [4] which considers that a-chains 
are the first to react. With this hypothesis, the rate constants of the different chemical reactions have been 
determined by multiple experimental measurements which references are available in [4|. The diagram 
of the reactions involved in the model is shown on figure [15] Note that because these experiments are not 
made at the temperature of human body, adjustments have been made by multiplication with a unique factor 
each of the constant rate of bi-molecular reactions. Once these adjustments done, we have successfully 
compared the resulting chemical model with known properties of hemoglobin. 
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Figure 15: Hemoglobin chemistry diagram (from [4]): hemoglobin has two allosteric states T and R (resp. 
exponent and 1). Hemoglobin can go from one state to the other and can carry from to 4 oxygen 
molecules. 

In the following, we will write M the concentration of hemoglobin of state j g {0, 1} (j = corre- 
sponds to T state while j — 1 corresponds to R state) carrying i e {0, 1, 2, 3, 4} oxygen molecules (noted 
Hi on figure [TBIl. 

An important number in blood studies is hemoglobin saturation, this number corresponds to the percent- 
age of hemoglobin capacity which is used to carry oxygen. For instance, during rest regime hemoglobin 
saturation is 0.97 after crossing the lungs while it is 0.75 before. It is the ratio between the quantity of O2 
effectively captured over the maximum quantity of O2 possible to capture. In the model chosen here, it is 
calculated by: 

1 + faf) + 2{h T 2 + fag) + S{hl + hj) + 4(hJ + faf) 

" U2 4 {h T Q + fa*) + {hi + faf ) + {h T 2 + fa?) + {hi + h§) + {hj + hi) 1 1 

On figure[l6j the saturation curve resulting from the model is shown. It has been plotted against oxygen 
partial pressure in the medium (mmHg). The steep part has an important role in oxygen exchange and 
its position is important, see [23] for more details. The model fits well measurements, for example the 
partial pressure in oxygen needed to have a 0.5 saturation (steep part of the curve) is 27 mmHg while it is 
measured to 26.8 mmHg in [23 j. 



3.2 Equations 

The path of oxygen from the aveolus into an hemoglobin molecule is described on figure [T7] In order to 
define the equations, we will note P the plasma region and G the cytosol. We will integrate membranes 
crossing in the boundary conditions in order to simplify the problem, we will note M the boundary between 
alveolus and capillary and F = dP\M. We will call v the plasma velocity which is considered as a data 
in this part, but which comes from calculations of the preceding sections. We will call a the oxygen partial 
concentration in plasma and a g the oxygen concentration in the cytosol G. 

All parameters used in the following equations comes from the references l4l [TTl[T7ll23l . 

Oxygen equations on P are convection/diffusion equations written in the deformed frame for all t > 0: 
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Figure 16: Hemoglobin saturation given by the model versus oxygen partial pressure in the medium 
(1 mmHg = 133 Pa). The curve fits well measured data given in the literature [23]. 
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Figure 17: Path of oxygen in a pulmonary capillary: oxygen first crosses the alveolus/capillary membrane 
dissolving into the plasma. Then it is transported in the plasma by convection/diffusion and reaches the red 
blood cell. It then crosses red blood cell membrane and diffuses inside the cytosol until it reacts with an 
hemoglobin molecule. 



|+I(«-^«)-V]o-fliAa = onP 
D 1 §±(Z,t) = ^(f3P A -a(Z,t))xl A (Z + m c (t).E c ) for Z E M and t > 
£) 1 |^ = onF {) 

1 an tg y " 1 

where D\ is the diffusion coefficient of oxygen in plasma, D2 the diffusion coefficient of oxygen in the 
tissues, r ac is the alveolus/capillary membrane thickness, (3 the solubility coefficient of oxygen in plasma 
and Pa the oxygen partial pressure in the aveolus. The supplementary term 1\(Z + m c (t).E c ) corresponds 
to the fact that the camera is moving and determines when it will reach a region permeable to oxygen. Thus, 
the set A corresponds to the portion of capillary which is permeable to oxygen. Typically, in our simulations 
: A = {X\500fj,m < X c < L c + 500/jto} with L c being the capillary typical length (1 mm), see section 
|4] We assume that no diffusive flow crosses the inlet and the outlet. Finally tq is the thickness of the red 
blood cell membrane and we recall that a g is the concentration of oxygen in the cytosol G. 

Note that because v is close to dm c /dt near the red blood cell, diffusion is the major phenomena 
involved in the transport of oxygen near the globule. However, close to the wall of capillary, the apparent 
transport velocity of oxygen is dm c /dt ~ 1 mm/ s (because v ~ 0) and this velocity is of the same order 
than the diffusive velocity D\/d c = 24.1CP 10 /2.10~ 6 ~ Imm/s (d c is the mean distance to be crossed 
by diffusion). Hence convection and diffusion each play a non negligible role in oxygen repartition near 
the capillary wall (the Peclet number is 1). 

We will now write hemoglobin and oxygen equations in the cytosol G. Each hemoglobin molecule Hi 
will diffuse in the cytosol (no transport) and can either: 
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• capture an oxygen molecule except if i = 4, at rate K\ . 

• free an oxygen molecule except if i = 0, at rate kf . 

• change its configuration, j — ► 1 — j, at rate c{. 

Hence, doing the balance of sources and sinks of each quantities, we can write, in the deformed frame, 
the following equations for H? concentration hj, for all j = 1.2 and t > 0: 

' °g - D 3 A hi = -Kja g hj 4h{ + ki +1 h{ +1 + c\ ■' on G, for i = 

< D 3 A^ -K{a g hi - hi h[ - 4h{ + Kl_ x a g hU + H+M +1 + c] V > on G, for i + 0, 4 

/■•//,; - «•;//; + KUa g h\_ x + ' on G, for, = 4 

, =° ondG,Vz = 0,...,4 

(9) 

Z3 3 is the hemoglobin diffusion coefficient in the cytosol. The boundary condition means that hemoglobin 
does not cross the red blood cell membrane. 

Oxygen diffuses into the cytosol and reacts with the different molecules of hemoglobin. The balance 
for oxygen in the cytosol can then be written in the deformed frame for all t > 0: 

^ lh - n ;! V. , ; A7«,/< • V . :: J/lr; onG 
D 2 ^ = ^(a-a g ) ondG 

We assume that the diffusion coefficient of oxygen in the cytosol is smaller than in plasma, because 
hemoglobin are big molecules. It has been estimated to be of the same order than in the tissue, see Ifl2l . 

The initial conditions are chosen in order to correspond to poor oxygenated blood with 40 ramHg of 
oxygen partial pressure. Initial hemoglobin state is assumed at equilibrium everywhere in the red blood 
cell with this oxygen partial pressure. 

4 Oxygen and hemoglobin: first numerical results 

We will restrict this preliminary study to two similar cases with pulmonary capillaries which geometries 
are simple bottlenecks as on figure [7] but with different sizes: the structure is 2 mm long and the con- 
stricted part, where oxygen exchanges will take place, is 1 mm long (typical length of a capillary see [23]). 
The larger diameter will be 8.5 fim in both cases. The smaller diameter will be 6.5 fim in the first case 
and 7.5 /j,m in the second one. These geometries will be crossed by five red blood cells. We assume that 
oxygen exchanges with alveolus stalls at the beginning of the bottleneck. 



On figure[l8j the mean saturation of oxygen is plotted against the mean position of the red blood cells 
in the capillary. The first case, with smaller diameter, induces a quicker capture of oxygen by hemoglobin. 
This is mainly due to the fact that: 

• the distance between the red blood cells and the capillary wall is smaller in the first case than in the 
second case, hence oxygen crosses the gap quicker when the diameter is small. 

• the red blood cell surface "available" for oxygen from capillary wall is larger when the diameter is 
smaller. Actually, for the larger diameter case, a large part of the surface stays in the middle of the 
capillary and oxygen needs to cross a larger distance to reach it. 

These observations are shown on figure [191 where the middle red blood cell is drawn for both cases. 
Both images are taken at the same abscissa, and their difference in term of oxygen saturation is important. 
As expected the red blood cell regions the more accessible to oxygen (closer to the capillary wall) are first 
saturated even though hemoglobin is diffusing inside the cytosol. Note that the less saturated region is on 
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Figure 18: Evolution of the mean hemoglobin saturation 5*02 in the red blood cells along the capillary. 
The capillary becomes permeable to oxygen at 500/jto. The larger the radius of the capillary, the slower 
hemoglobin captures oxygen. Note that because of plasma flow effect on oxygen distribution in the capil- 
lary, some slight SO2 inhomogeneities arises before the capillary becomes permeable. 



the axis, in the rear part of the red blood cell where oxygen flow is the smaller. This effect is due to the 
screening made in this region by the red blood cell shape to the diffusive process. The screening effect of 
diffusion can thus play an important role in biological systems. This influence has already been shown in 
the case of the lung, see (8). Note that the presence of the other four red blood cells also limits oxygen 
capture speed through the same phenomena. 
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Figure 19: Saturation of oxygen in the middle red blood cell. The arrows represent the flow of oxygen in 
plasma. The upper image corresponds to a 6.5 ixm diameter and the lower image corresponds to a 7.5 fim 
diameter. These images are taken near a 1300/im progression in the capillary. Note the differences in 
oxygen saturation of hemoglobin. 



Conclusion 

In this paper, a full modeling chain of red blood cell and oxygen behaviors in a capillary has been built. This 
gives access to a large panel of data which should be useful in order to better understand oxygen transfer 
at capillary level. For instance, the fine tuning of some physiological parameters (for example as Vo 2 max 
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in the lungs or as blood speed regulation in the vascular system) should be explained by the phenomena 
involved in capillaries and help to understand the evolution processes that led to such structures. Moreover, 
pathological behavior (for example paludism or diabetes) could be integrated in this model and could help 
to understand the consequences of such dysfunctions. 

The prospects of this approach is to develop a three dimensional model and to add collisions between 
red blood cells and with capillary walls. There also exists an important interaction between erythrocytes: 
they can combine together to form rolls, up to eight cells. These rolls are maintained by chemical forces 
which can break under high shears. The integration of this behavior is important because of its role in blood 
circulation. Moreover the restriction to straight capillary could be removed. Finally, mathematical studies 
of the method should be performed in order to control precisely the errors made during the modeling. 
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